Efficient Square-Root Free 2 Symbol Max-Log Receiver

ABSTRACT

The present invention employs a look up table based implementation for the metric computations which eliminate redundancy and substantially reduce the number of multiplications. Moreover, inventive method exploits the fact that the un-normalized constellation symbols are complex integers so that the product of a real-number and an un-normalized constellation symbol can be implemented by only additions. The inventive method also enables a greater efficiency for whitening colored noise prior to demodulation, one of which involves no square-root operation. The invention results in less complexity, faster operation, lower power consumption, without sacrificing performance

RELATED APPLICATION INFORMATION

This application claims priority to provisional application serial number 61,469,240 filed on Mar. 30, 2011, the contents thereof are incorporated herein by reference.

BACKGROUND

The present invention relates to wireless communications and, more particularly, to an efficient square-root free 2 symbol max-log receiver.

In several communication scenarios the following receiver (demodulator) design problem is ubiquitous. Consider a received signal vector z that has 2 complex-valued elements and can be expressed as:

z=Hs+n  (Equation)

where, s is the transmit symbol vector having 2 complex-valued elements, each drawn from a normalized constellation. The matrix H models the channel and the vector n models the additive independent noise, assumed to have complex Gaussian elements of unit variance. Then given z and H and the constellation(s) from which the elements of s are drawn, the receiver (demodulator) design problem is to determine the optimal hard decision about s and/or the optimal soft decisions (log-likelihood ratios) about the coded bits mapped to s.

One approach to solving the above problem, the brute-force way, to determine either the hard or soft decisions is to list out all possibilities of s and compute associated metrics. This method has a very high complexity which scales as O(M²), where M is the constellation size and is considered to be impractical. A better approach was consequently developed in U.S. patent application Ser. No. 11/857,269, inventors: Prasad et al., entitled “Max_Log Receiver For Multiple Input Multiple output (MIMO) Systems”. In this technique, the demodulator involves twice linearly transforming the received signal vector to obtain new transformed vectors that can be modeled as in the above Equation, but where the transformed channel matrices have a triangular structure. The transformed vectors are then used for metric computations after exploiting the induced triangular structure in the transformed channel matrices. Determining the matrices used for these two linear transformations as well as the elements of the resulting transformed channel matrices involves square-root operations that are costly to implement.

Accordingly, there is a need for an efficient square-root free 2 symbol max-log receiver.

SUMMARY

A method for a square-root free 2 symbol max-log receiver includes obtaining linear transformations of a received two stream signal and a channel matrix without implementing square-root operations, listing out all possibilities for a first symbol of the received two stream signal, building look-up-tables needed for computing first metrics associated with all possibilities for a first symbol of the two stream signal, determining a second symbol of the two stream signal for each the first symbol listed out, evaluating said first metrics for each the first symbol and second symbol pair using the look up tables, listing out all possibilities for the second symbol 2, building look-up-tables needed for computing second metrics associated with all possibilities for a second symbol of said received two stream signal, determining a first symbol for each choice of the second symbol listed out, evaluating the second metrics for each the second symbol and first symbol pair using the look up tables, determining an exact max-log log likelihood ratio for each coded bit using the second metrics; and decoding a least one codeword in the two stream signal using the determined exact max-log log likelihood ratios for all bits.

These and other features and advantages will become apparent from the following detailed description of illustrative embodiments thereof, which is to be read in connection with the accompanying drawings.

BRIEF DESCRIPTION OF DRAWINGS

The disclosure will provide details in the following description of preferred embodiments with reference to the following figures wherein:

FIG. 1 is a diagram of an exemplary system schematic in which the present inventive principles can be practiced.

FIG. 2 is a block diagram of a receiver configuration in accordance with the invention.

FIG. 3 is a flow diagram for an efficient square-root free 2 symbol max-log receiver, in accordance with the invention.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

Referring now in detail to the figures in which like numerals represent the same or similar elements and initially to FIG. 1, is a diagram of an exemplary system schematic in which the present inventive principles can be practiced. A destination receiver 10 receives signals of interest 11, 12 as well as interfering signals. This invention improves upon the prior method noted above by employing a new way of linear transformations and metric computations that are square-root free.

The present invention employs a look up table based implementation for the metric computations which eliminate redundancy and substantially reduce the number of multiplications. Moreover, inventive method exploits the fact that the un-normalized constellation symbols are complex integers so that the product of a real-number and an un-normalized constellation symbol can be implemented by only additions. The inventive method also enables a greater efficiency for whitening colored noise prior to demodulation, one of which involves no square-root operation. The invention results in less complexity, faster operation, lower power consumption, without sacrificing performance

Referring now to FIG. 2, the block diagram shows a receiver configuration in accordance with the invention. The block 20 coupled to receive antennas 1 . . . Nr provides base-band converted output to the inventive efficient modulator 21 that is influenced by a channel and interference covariance matrix estimator 22. The key ingredients of the present invention are square-root free operation and look-up table based implementation.

Detailed Process

Signal Model:

$z = {{{H\begin{bmatrix} E^{(0)} & 0 \\ 0 & E^{(1)} \end{bmatrix}}\begin{bmatrix} x^{(0)} \\ x^{(1)} \end{bmatrix}} + \begin{bmatrix} n^{(0)} \\ n^{(1)} \end{bmatrix}}$ $z = {{{{H\begin{bmatrix} E^{(0)} & 0 \\ 0 & E^{(1)} \end{bmatrix}}x} + n} = {{\overset{\_}{H}x} + n}}$

Equation 1: Rx Signal Model for Spatial Multiplexing

Where H is effective channel matrix of size N_(r)x2, N_(r) being the number of receive antennas, n˜N_(c)(0,I) and x is a vector of un-normalized QAM symbols. E⁽⁰⁾, E⁽¹⁾ are normalizing factors such that E⁽⁰⁾x⁽⁰⁾, E⁽¹⁾x⁽¹⁾ have unit average energy. Thus, x⁽⁰⁾, x⁽¹⁾ are un-normalized QAM symbols.

Given

Received vector z,

Estimate {tilde over (H)}=[{tilde over (h)}₀, {tilde over (h)}₁] of effective channel matrix H

Using QR decomposition, {tilde over (H)} can be factored as

{tilde over (H)}=QR  Equation 2

Where Q is N_(r)x2 semi-unitary matrix i.e. Q^(H)Q=I and R is 2×2 upper triangular matrix that can be expanded as

$R = \begin{bmatrix} r_{0,0} & r_{0,1} \\ 0 & r_{1,1} \end{bmatrix}$

Square-Root Free Operation

Referring now to the flow diagram of FIG. 3 showing an efficient square-root free 2 symbol max-log receiver, in accordance with the invention, the square-root free linear transformations of the received signal and the channel matrix 100 is obtained as described below.

Obtain

${\overset{\sim}{Q} = {Q\begin{bmatrix} {1/r_{0,0}} & 0 \\ 0 & {1/r_{1,1}} \end{bmatrix}}},{{\overset{\Cup}{r}}_{0,1} = {r_{0,1}/r_{0,0}}}$ and θ = r_(1, 1)²/r_(0, 0)²

Upon multiplying received vector z with {tilde over (Q)}^(H) we have

$\begin{matrix} {\mspace{79mu} {{\overset{\sim}{z} = {{{\overset{\sim}{Q}}^{H}z} = {{{\overset{\sim}{Q}}^{H}\left( {{QRx} + n} \right)} = {{\overset{\sim}{R}x} + \overset{\sim}{n}}}}}{\overset{\sim}{z} = {\begin{bmatrix} {\overset{\sim}{z}}^{(0)} \\ {\overset{\sim}{z}}^{(1)} \end{bmatrix} = {{{\begin{bmatrix} 1 & {r_{0,1}/r_{0,0}} \\ 0 & 1 \end{bmatrix}\begin{bmatrix} x^{(0)} \\ x^{(1)} \end{bmatrix}} + \overset{\sim}{n}} = {\begin{bmatrix} {x^{(0)} + {{\overset{\Cup}{r}}_{0,1}x^{(1)}}} \\ x^{(1)} \end{bmatrix} + \begin{bmatrix} {\overset{\sim}{n}}^{(0)} \\ {\overset{\sim}{n}}^{(1)} \end{bmatrix}}}}}}} & {{Equation}\mspace{14mu} 3} \end{matrix}$

Note that since Q^(H)Q=I=>ñ˜N_(c)(0, diag{1/r_(0,0) ², 1/r_(1,1) ²}).

To implement the inventive procedure, the term {tilde over (Q)}, {hacek over (r)}_(0,1), r_(1,1) ² and θ needs to be determined. Also note that r_(1,1) ² and θ are both real and positive valued, whereas {hacek over (r)}_(0,1) is complex valued. The computation of these quantities is described in the following. Note that no square-root operations are involved.

Obtaining {tilde over (Q)},{hacek over (r)}_(0,1), r_(1,1) ² and θ via a modified QR decomposition (For any N_(r)≧2), use {tilde over (H)} in the following steps

-   -   Step 1: Consider {tilde over (H)}=[{tilde over (h)}₀, {tilde         over (h)}₁] and obtain

α={tilde over (h)} ₀ ^(H) {tilde over (h)} ₁

-   -   Step 2: Obtain δ=∥{tilde over (h)}₀∥² and first column of {tilde         over (Q)}=[{tilde over (q)}₀, {tilde over (q)}₁] as

{tilde over (q)} ₀ ={tilde over (h)} ₀/δ

-   -   Step 3: Obtain the vector {tilde over (v)}={tilde over         (h)}₁−α{tilde over (q)}₀; and r_(1,1) ²=∥{tilde over (v)}∥² and         the second column of {tilde over (Q)}=[{tilde over (q)}₀,{tilde         over (q)}₁] as

{tilde over (q)} ₁ ={tilde over (v)}/r _(1,1) ²

-   -   Step 4: Obtain {hacek over (r)}_(0,1) as

{hacek over (r)} _(0,1)=α/δ

-   -   Step 5: Obtain θ as

θ=r _(1,1) ²/δ

Referring now to FIG. 3 and blocks 103 and 106, the following are denoted:

-   -   K^((v)) as number of symbols in modulation constellation applied         to layer v. K^((v)) has value 4,16 and 64 if modulation is QPSK,         16QAM and 64QAM respectively     -   S^((v)) as the set containing K^((v)) modulation symbols in the         modulation constellation applied to layer v.     -   B^((v)) as a number of bits per modulation symbol applied to         layer v. B^((v)) has value 2,4 and 6 if modulation is QPSK,         16QAM and 64QAM respectively

The soft bit k:k=0, 1, . . . B^((v))−1, correspond to transmitted symbol x^((v)): v=0, 1 are generated using max-log approximation of LLR as follows:

$\begin{matrix} {{L\left( b_{k,v} \right)} \approx {{\min\limits_{c_{j} \in S_{k,v,0}}{{\left( {\overset{\sim}{z} - {\overset{\sim}{R}c_{j}}} \right)^{H}\begin{bmatrix} r_{0,0}^{2} & 0 \\ 0 & r_{1,1}^{2} \end{bmatrix}}\left( {\overset{\sim}{z} - {\overset{\sim}{R}c_{j}}} \right)}} - {\min\limits_{c_{j} \in S_{k,v,1}}{{\left( {\overset{\sim}{z} - {\overset{\sim}{R}c_{j}}} \right)^{H}\begin{bmatrix} r_{0,0}^{2} & 0 \\ 0 & r_{1,1}^{2} \end{bmatrix}}\left( {\overset{\sim}{z} - {\overset{\sim}{R}c_{j}}} \right)}}}} & {{Equation}\mspace{14mu} 4} \end{matrix}$

Where

-   -   Vector c_(j)=[c_(j) ⁽⁰⁾, c_(j) ⁽¹⁾]^(T) with each element c_(j)         ^((v)) being a complex symbol in the modulation constellation         applied to layer v.     -   S_(k,v,1) and S_(k,v,0) denote the set containing all possible         vector c_(j) such that bit k of element c_(j) ^((v)) is equal to         1 and 0 respectively.         For example, if QPSK is used on layer v=0 and 16QAM is used on         layer v=1 then     -   S_(k,0,1) and S_(k,0,0) each contains 2×16=32 vectors     -   S_(k,1,1) and S_(k,1,0) each contains 8×4=32 vectors

Since full computation according to Equation 4 results in undesirably high complexity, further simplification is needed and is explained below. Note that the following simplification results in no loss of optimality when compared to (Equation 4). In other words, the soft bits generated below are identical respectively to the ones generated using (Equation 4).

Considering soft bit calculation for the second layer (i.e. v=1), following steps are performed

Step1: for each c_(j) ⁽¹⁾εS⁽¹⁾:j=0, 1, . . . , K⁽¹⁾−1,

-   -   Compute d_(j) ⁽¹⁾=θ|{tilde over (z)}⁽¹⁾−c_(j) ⁽¹⁾|⁽²⁾ or (see         Equation 3)     -   Compute soft estimate {tilde over (z)}⁽⁰⁾−{hacek over         (r)}_(0,1)c_(j) ⁽¹⁾ (see Equation 3) and then select a         symbol)ĉ_(j) ⁽⁰⁾εS⁽⁰⁾ that is closest to {tilde over         (z)}⁽⁰⁾−{hacek over (r)}_(0,1)c_(j) ⁽¹⁾     -   Compute d_(j) ⁽⁰⁾=|{tilde over (z)}⁽⁰⁾−ĉ_(j) ⁽⁰⁾−{hacek over         (r)}_(0,1)c_(j) ^((1)|) ²     -   Compute total distance d_(j) ^(TOTAL)=d_(j) ⁽⁰⁾+d_(j) ⁽¹⁾

Step2:

-   -   Determine

${{\hat{d}}^{(1)} = {\min\limits_{c_{j}^{(1)} \in S^{(1)}}\left\{ d_{j}^{TOTAL} \right\}}},{{\hat{c}}_{j}^{(1)} = {\underset{c_{j}^{(1)} \in S^{(1)}}{\arg \; \min}\left\{ d_{j}^{TOTAL} \right\}}}$

and note that ĉ_(j) ⁽¹⁾εS⁽¹⁾ is the symbol of layer−1 in the maximum likelihood (ML) decision. Let {circumflex over (b)}_(k,1): k=0, 1, . . . , B⁽¹⁾−1 be the bit-label of ĉ_(j) ⁽¹⁾εS⁽¹⁾

-   -   Compute each soft bit L(b_(k,1)): k=0, 1, . . . , B⁽¹⁾−1 as

$\begin{matrix} {{{L\left( b_{k,1} \right)} = {{{r_{0,0}^{2}\left( {{\hat{d}}^{(1)} - {\min\limits_{d_{m}^{TOTAL} \in D_{k}^{1}}d_{m}^{TOTAL}}} \right)}\mspace{14mu} {if}\mspace{14mu} {\hat{b}}_{k,1}} = 0}}{or}{{{L\left( b_{k,1} \right)} = {{{r_{0,0}^{2}\left( {{\min\limits_{d_{m}^{TOTAL} \in D_{k}^{0}}d_{m}^{TOTAL}} - {\hat{d}}^{(1)}} \right)}\mspace{14mu} {if}\mspace{14mu} {\hat{b}}_{k,1}} = 1}};}} & {{Equation}\mspace{14mu} 5b} \end{matrix}$

-   -   Where D_(k) ¹ and D_(k) ⁰ denotes the set containing all         distances d_(j) ^(TOTAL) correspond to c_(j) ⁽¹⁾ that has that         bit k equal to 1 and 0 respectively.

Considering soft bit calculation for the first layer (i.e. v=0), following steps are performed

-   -   Step A: Swapping of symbol vector x and estimate of effective         channel matrix {tilde over (H)} as

$x^{\prime} = \begin{bmatrix} x^{(1)} \\ x^{(0)} \end{bmatrix}$ and ${\overset{\sim}{H}}^{\prime} = \left\lbrack {{\overset{\sim}{h}}_{1},{\overset{\sim}{h}}_{0}} \right\rbrack$

-   -   Step B: Perform QR decomposition as shown in Equation 2 (using         {tilde over (H)}′ instead of {tilde over (H)}) and multiply         receive vector with {tilde over (Q)}^(H) as shown in Equation 3     -   Step1 and Step2 are similar to description above with         appropriate swapping i.e. S⁽¹⁾, S⁽⁰⁾, K⁽¹⁾, B⁽¹⁾, L(b_(k,1))         being replaced by S⁽⁰⁾, S⁽¹⁾, K⁽⁰⁾, B⁽⁰⁾, L(b_(k,0))         respectively.

Attention is now directed to Look-Up-Table (LUT) based metric computations denoted in blocks 102 and 105 of FIG. 3. In this section we simplify the computations involved in determining d_(j) ⁽⁰⁾, d_(j) ⁽¹⁾. This is particularly important when at-least one of the two constellations is 64 QAM. The key ideas we use are the following:

-   -   1. We look at the max-log LLR expression and after expanding         d_(j) ⁽⁰⁾, d_(j) ⁽¹⁾, we drop the terms that do not influence         the LLR.     -   2. In the model in (Equation 1) we pull out normalizing factors         so that the modulation symbols are from un-normalized QAM         constellations which contain integers. Then we can use the fact         that the product of an integer and a real number can be obtained         using just additions.     -   3. We also form LUTs containing common terms that mostly depend         only on the channel coefficients. The entries of these LUTs are         repeatedly accessed while computing all d_(j) ⁽⁰⁾, d_(j) ⁽¹⁾ and         multiplications are avoided.

Consider the step: Compute d_(j) ⁽¹⁾=θ|{tilde over (z)}⁽¹⁾−c_(j) ⁽¹⁾|²

Note that we can expand it as

d _(j) ⁽¹⁾ =θ|{tilde over (z)} ⁽¹⁾|²+θ(c _(j,R) ⁽¹⁾)²+θ(c _(j,I) ⁽¹⁾)²−2{tilde over (z)} _(R) ⁽¹⁾θ−2{tilde over (z)} _(I) ⁽¹⁾ c _(j,I) ⁽¹⁾θ  (Equation 10)

From (Equation 4b), we can conclude that computing (Equation 10) instead as

d _(j) ⁽¹⁾=θ(c _(j,R) ⁽¹⁾)²+θ(c _(j,I) ⁽¹⁾)²−2{tilde over (z)} _(R) ⁽¹⁾ c _(j,R) ⁽¹⁾θ−2{tilde over (z)} _(I) ⁽¹⁾ c _(j,I) ⁽¹⁾θ  (Equation 11)

results in no loss of max-log LLR optimality. Then since c_(j) ⁽¹⁾ belongs to the un-normalized QAM constellation S⁽¹⁾=S_(R) ⁽¹⁾+iS_(I) ⁽¹⁾, where S_(R) ⁽¹⁾, S_(I) ⁽¹⁾ are both identical un-normalized PAM constellations of size M⁽¹⁾=√{square root over (K⁽¹⁾)} and given by {−(M⁽¹⁾−1), −(M⁽¹⁾−3), . . . , −1, 1, . . . , (M⁽¹⁾−3), (M⁽¹⁾−1)}. Then in a one dimensional (1-D) LUT of length M⁽¹⁾/2 we can pre-compute and store {(M⁽¹⁾−1)²θ, (M⁽¹⁾−3)²θ, . . . , θ}. Note that since the product of any positive integer and any real number can be determined only by additions, this 1-D LUT can be determined only by additions and needs to be updated only if θ changes. Each c_(j,R) ⁽¹⁾ and c_(j,I) ⁽¹⁾ should then index the appropriate entry of this LUT via a pre-defined mapping rule to obtain (c_(j,R) ⁽¹⁾)²θ and (c_(j,I) ⁽¹⁾)²θ, respectively. Further, two other 1-D LUTs containing {2(M⁽¹⁾−1){tilde over (z)}_(R) ⁽¹⁾θ,2(M⁽¹⁾−3){tilde over (z)}_(R) ⁽¹⁾θ, . . . , 2{tilde over (z)}_(R) ⁽¹⁾θ} and {2(M⁽¹⁾−1){tilde over (z)}_(I) ⁽¹⁾θ,2(M⁽¹⁾−3){tilde over (z)}_(I) ⁽¹⁾θ, . . . , 2{tilde over (z)}_(I) ⁽¹⁾θ}, respectively, can be generated via only addition operations, once the two terms {tilde over (z)}_(R) ⁽¹⁾θ,{tilde over (z)}_(I) ⁽¹⁾θ are computed. Then all possibilities of the terms −2{tilde over (z)}_(R) ⁽¹⁾c_(j,R) ⁽¹⁾θ−2{tilde over (z)}_(I) ⁽¹⁾c_(j,I) ⁽¹⁾θ can be determined using these two 1-D LUTs by accessing the proper entries followed by a negation if needed. Alternatively, once the two terms {tilde over (z)}_(R) ⁽¹⁾θ,{tilde over (z)}_(I) ⁽¹⁾θ are computed, each −2{tilde over (z)}_(R) ⁽¹⁾c_(j,R) ⁽¹⁾θ−2{tilde over (z)}_(I) ⁽¹⁾c_(j,I) ⁽¹⁾θ can be directly computed using additions. As a result, all {d_(j) ⁽¹⁾} can be computed with only 2 real multiplications.

Next, consider the computation of

d _(j) ⁽⁰⁾ =|{tilde over (z)} ⁽⁰⁾ −ĉ _(j) ⁽⁰⁾ −{hacek over (r)} _(0,1) c _(j) ⁽¹⁾|²=({tilde over (z)} _(R) ⁽⁰⁾ −ĉ _(j,R) ⁽⁰⁾ −{hacek over (r)} _(0,1,R) c _(j,R) ⁽¹⁾ +{hacek over (r)} _(0,1,I) c _(j,I) ⁽¹⁾)²+({tilde over (z)} _(I) ⁽⁰⁾ −ĉ _(j,I) ⁽⁰⁾ −{hacek over (r)} _(0,1,I) c _(j,R) ⁽¹⁾ −{hacek over (r)} _(0,1,R) c _(j,I) ⁽¹⁾)²  (Equation 12)

Consider the first term in (Equation 12) which is ({tilde over (z)}_(R) ⁽⁰⁾−ĉ_(j,R) ⁽⁰⁾−{hacek over (r)}_(0,1,R)c_(j,R) ⁽¹⁾+{hacek over (r)}_(0,1,I)c_(j,I) ⁽¹⁾)². Suppose that we have computed via simple quantizing the terms {{tilde over (z)}_(R) ⁽⁰⁾−ĉ_(j,R) ⁽⁰⁾−{hacek over (r)}_(0,1,R)c_(j,R) ⁽¹⁾+{hacek over (r)}_(0,1,I)c_(j,I) ⁽¹⁾} for all c_(j) ⁽¹⁾=c_(j,R) ⁽¹⁾+ic_(j,I) ⁽¹⁾. Note that there are K⁽¹⁾ such terms and our objective is to compute the squares of these terms as efficiently as possible. The obvious method would entail K⁽¹⁾ real multiplications. Here we propose an LUT based efficient method and to illustrate, we consider 64 QAM so that c_(j,R) ⁽¹⁾ and c_(j,I) ⁽¹⁾ each belong to the set {−7, −5, . . . , −1, 1, . . . , 5, 7}. Then, suppose that for any fixed c_(j) ₁ ⁽¹⁾=a_(R)−i7, we have computed the term w=γ_(R,j)−ĉ_(j) ₁ _(,R) ⁽⁰⁾−7{hacek over (r)}_(0,1,I) and w², where γ_(R,j)={tilde over (z)}_(R) ⁽⁰⁾−{hacek over (r)}_(0,1,R)a_(R). Moreover, suppose that for any c_(j) ₂ ⁽¹⁾=a_(R)+iq, q={−5, −3, . . . , 7} we have determined u=γ_(R,j)−ĉ_(j) ₂ _(,R) ⁽⁰⁾+q{hacek over (r)}_(0,1,I). Then letting {circumflex over (d)}=ĉ_(j) ₂ _(,R) ⁽⁰⁾−ĉ_(j) ₁ _(,R) ⁽⁰⁾ we see that u² can be expanded as

u ² =w ²+2(7+q)w{hacek over (r)} _(0,1,I)−2{circumflex over (d)}w+{circumflex over (d)} ²+((7+q){hacek over (r)} _(0,1,I))²−2{circumflex over (d)}(7+q){hacek over (r)} _(0,1,I)  (Equation 13)

Note in (Equation 13) that if w², {hacek over (r)}_(0,1,I), ({hacek over (r)}_(0,1,I))², w{hacek over (r)}_(0,1,I) are available then u² can be determined via only addition (or subtraction which can be considered equivalent to an addition) operations since {circumflex over (d)} is an integer and the product of an integer, say 3, and a real number, say x, is equal to x+x+x.

To efficiently determine any such u², we construct a 2-D LUT which for a given pair {circumflex over (d)}, q yields {circumflex over (d)}²+((7+q){hacek over (r)}_(0,1,I))²−2{circumflex over (d)}(7+q){hacek over (r)}_(0,1,I). Since all possible {circumflex over (d)}, q are integer pairs, this 2-D LUT can be constructed using only addition operations and needs to be updated only if {hacek over (r)}_(0,1,I) changes. Furthermore, for a given pair w, w{hacek over (r)}_(0,1,I) the quantities {2(7+q)w{hacek over (r)}_(0,1,I)} and {2{circumflex over (d)}w} can either be directly computed using additions or can be recursively updated (again using additions) as q is varied from −5 to 7 in steps of 2. Thus, for a given {hacek over (r)}_(0,1,I), ({hacek over (r)}_(0,1,I))² only two real multiplications (i.e., those needed to compute w², w{hacek over (r)}_(0,1,I)) are required to determine all {({tilde over (z)}_(R) ⁽⁰⁾−ĉ_(j,R) ⁽⁰⁾−{hacek over (r)}_(0,1,R)c_(j,R) ⁽¹⁾+{hacek over (r)}_(0,1,I)c_(j,I) ⁽¹⁾)²}, c_(j,R) ⁽¹⁾=a_(R), c_(j,R) ⁽¹⁾=q for any fixed a_(R) as q is varied from −5 to 7 in steps of 2. It can be similarly shown that given {hacek over (r)}_(0,1,R), ({hacek over (r)}_(0,1,R))² only one more additional real multiplication (needed to compute w{hacek over (r)}_(0,1,R)) is required to determine all {({tilde over (z)}_(R) ⁽⁰⁾−ĉ_(j,R) ⁽⁰⁾−{hacek over (r)}_(0,1,R)c_(j,R) ⁽¹⁾+{hacek over (r)}_(0,1,I)c_(j,I) ⁽¹⁾)²}, c_(j,R) ⁽¹⁾=a_(R), c_(j,I) ⁽¹⁾=q for any fixed a_(R) as q is varied from −5 to 7 in steps of 2. Thus, all {({tilde over (z)}_(R) ⁽⁰⁾−ĉ_(j,R) ⁽⁰⁾−{hacek over (r)}_(0,1,R)c_(j,R) ⁽¹⁾+{hacek over (r)}_(0,1,I)c_(j,I) ⁽¹⁾)²} for all c_(j) ⁽¹⁾ can be determined with 10 real multiplications. Similarly all possibilities of the second term in (Equation 12) which are {({tilde over (z)}_(I) ⁽⁰⁾−ĉ_(j,I) ⁽⁰⁾−{hacek over (r)}_(0,1,I)c_(j,R) ⁽¹⁾+{hacek over (r)}_(0,1,R)c_(j,I) ⁽¹⁾)²} be determined with 10 real multiplications. Consequently all {d₀ ⁽¹⁾} can be computed with only 20 real multiplications.

Detailed Implementation (for 2 RX)

1.1 Second Layer Processing

Set

-   -   v=0     -   S^(A)=S⁽¹⁾     -   S^(B)=S⁽⁰⁾     -   K^(A)=K⁽¹⁾     -   B^(A)=B⁽¹⁾         Then perform processing as described in subsections below.         After completing processing in 1.1.4, assign L(b_(k,1))=L(b_(k))

1.1.1 Process Received Input

$\begin{matrix} {\overset{\sim}{z} = {\begin{bmatrix} {\overset{\sim}{z}}^{(0)} \\ {\overset{\sim}{z}}^{(1)} \end{bmatrix} = {{\overset{\sim}{Q}}^{H}z}}} & {{Equation}\mspace{14mu} 6} \end{matrix}$

Obtain {hacek over (r)}_(0,1), r_(1,1) ² and θ.

1.1.2 Distance Calculation 1

For each c_(j) ^(A)εS^(A):j=0, 1, . . . , K^(A)−1,

Step 1: Compute

d _(j) ^(A) =θ|{tilde over (z)} ⁽¹⁾ −c _(j) ^(A)|²  Equation 7

Note that the computation in Equation 15 can be simplified using LUT as described above

Step 2: Compute the soft estimate

{tilde over (c)} _(j) ^(B) ={tilde over (z)} ⁽⁰⁾ −{hacek over (r)} _(0,1) c _(j) ^(A)  Equation 8

and compute ĉ_(j) ^(B) which is the symbol in S^(B) closest to a {tilde over (c)}_(j) ^(B) via quantization.

1.1.2 Distance Calculation 2

Step 1: Compute

d _(j) ^(B) =|{hacek over (z)} ⁽⁰⁾ −ĉ _(j) ^(B) −{hacek over (r)} _(0,1) c _(j) ^(A)|² :j=0, 1, . . . , K ^(A)−1  Equation 9

Note that the computation in Equation 17 can be significantly simplified using LUTs as described above

Step 2: Compute

d _(j) ^(TOTAL) =d _(j) ^(B) +d _(j) ^(A) :j=0, 1, . . . , K ^(A)−1  Equation 10

-   -   Determine

${{\hat{d}}^{TOTAL} = {\min\limits_{{j = 0},1,\mspace{11mu} \ldots \mspace{14mu},{K^{A} - 1}}\left\{ d_{j}^{TOTAL} \right\}}},{\hat{j} = {\underset{{j = 0},1,\mspace{11mu} \ldots \mspace{14mu},{K^{A} - 1}}{\arg \; \min}\left\{ d_{j}^{TOTAL} \right\}}}$

and let {circumflex over (b)}_(k):k=0, 1, . . . , B^(A)−1 be the bit-label of c_(ĵ) ^((A)).

Soft Bit Calculation

Compute each soft bit L(b_(k)): k=0, 1, . . . , B^(A)−1 as

$\begin{matrix} {{{L\left( b_{k} \right)} = {{{r_{0,0}^{2}\left( {{\hat{d}}^{TOTAL} - {\min\limits_{d_{m}^{TOTAL} \in D_{k}^{1}}d_{m}^{TOTAL}}} \right)}\mspace{14mu} {if}\mspace{14mu} {\hat{b}}_{k}} = 0}}{or}{{L\left( b_{k} \right)} = {{{r_{0,0}^{2}\left( {{\min\limits_{d_{m}^{TOTAL} \in D_{k}^{0}}d_{m}^{TOTAL}} - {\hat{d}}^{TOTAL}} \right)}\mspace{14mu} {if}\mspace{14mu} {\hat{b}}_{k}} = 1}}} & {{Equation}\mspace{14mu} 11} \end{matrix}$

Where D_(k) ¹ and D_(k) ⁰ denotes the set containing all distances d_(j) ^(TOTAL) correspond to c_(j) ^(A) that has bit k equal to 1 and 0 respectively.

First Layer Processing

There is pre-processing step which involve swapping of symbol vector x and estimate of effective channel matrix {tilde over (H)} to obtain as vector x′ and matrix {tilde over (H)}′ as follow

$\begin{matrix} {{{x^{\prime} = \begin{bmatrix} x^{(1)} \\ x^{(0)} \end{bmatrix}};}{{\overset{\sim}{H}}^{\prime} = \begin{bmatrix} {\overset{\sim}{h}}_{1,1} & {\overset{\sim}{h}}_{1,0} \\ {\overset{\sim}{h}}_{0,1} & {\overset{\sim}{h}}_{0,0} \end{bmatrix}}} & {{Equation}\mspace{14mu} 20} \end{matrix}$

Set

-   -   v=1     -   S^(A)=S⁽⁰⁾     -   S^(B)=S⁽¹⁾     -   K^(A)=K⁽⁰⁾     -   B^(A)=B⁽⁰⁾         Then perform processing as described in subsections below.         After completing processing in 1.2.4, assign         L(b_(k,0))=L(b_(k)).

Computations Via QR Decomposition

As described before with {tilde over (H)} being replaced by {tilde over (H)}′. Note that some computations made for deriving the first QR decomposition can be re-used. Process received input: As in 1.1.1 with x being replaced by x′.

Distance Calculation 1: As in 1.1.2.

Distance calculation 2: As in [0037]. Soft bit calculation: As in 1.1.4

Noise-whitening

Consider the received signal model at UE of interest having index 0 and let UE with index 1 be the co-scheduled user.

$\begin{matrix} {y_{0} = {{H_{0}V_{0}x_{0}} + {H_{0}V_{1}x_{1}} + n_{0}}} \\ {{= {{{{\overset{\sim}{H}}_{0}x_{0}} + {{\hat{H}}_{0}x_{1}} + {n_{0}y_{0}}} \in C^{2 \times 1}}},{n_{0} \in C^{2 \times 1}},} \\ {{{H_{0} \in C^{2 \times 4}},{V_{k} \in C^{4 \times 2}},{x_{k} \in C^{2 \times 1}}}} \\ {{{{\overset{\sim}{H}}_{0} \in C^{2 \times 2}},{{\hat{H}}_{0} \in C^{2 \times 2}}}} \end{matrix}$

Assuming that UE-0 can estimate {tilde over (H)}₀, Ĥ₀, we will derive the whitening process. The whitening process can be done using a Cholesky decomposition of the noise-plus-interference covariance matrix.

In particular, we compute the covariance matrix as C=E[(Ĥ₀x₁+n₀)(Ĥ₀x₁+n₀)*]=I+Ĥ₀Ĥ₀* where we have assumed E[x₁x₁*]=E[n₀n₀*]=I. Note that if we have estimated Ĥ₀ we can use it to compute C or we can determine C using covariance estimation. We can then determine the Cholesky decomposition of C as C=LL* where L is a lower or upper triangular matrix with positive diagonal elements. We then determine L⁻¹ and L⁻¹y₀ which can be expanded as z₀=L⁻¹y₀=L⁻¹{tilde over (H)}₀x₀+ñ and note that E[ññ*]=I. We can use the two-symbol demodulator derived above on z₀ with L⁻¹{tilde over (H)}₀ being the effective channel matrix.

In the special case when C is a 2×2 matrix we can obtain the following closed form expressions: Let

$C = \begin{bmatrix} a & c \\ b & d \end{bmatrix}$

and note that since C is positive definite we must have a>0, d>0 and b=c*. Let L be lower triangular with positive diagonal elements. Then we can determine elements of L as

${L = \begin{bmatrix} l_{11} & 0 \\ l_{21} & l_{22} \end{bmatrix}},{l_{11} = \sqrt{a}},{l_{21} = {b/l_{11}}},{l_{22} = {\sqrt{d - {l_{21}}^{2}}.}}$

The inverse of L can now be computed as

$L^{- 1} = {\begin{bmatrix} {1/l_{11}} & 0 \\ {{- l_{21}}/\left( {l_{11}l_{22}} \right)} & {1/l_{22}} \end{bmatrix}.}$

Square-Root Free Noise-Whitening

Consider again the received signal model:

y ₀ ={tilde over (H)} ₀ x ₀ +ñ ₀

E[ñ ₀ ñ ₀ *]=C  (Equation 21)

and suppose

$C = \begin{bmatrix} a & c \\ b & d \end{bmatrix}$

is a 2×2 positive definite matrix. Note that

$C^{- 1} = {{\frac{1}{{ad} - {bc}}\begin{bmatrix} d & {- c} \\ {- b} & a \end{bmatrix}}.}$

Then obtain the following variation of the Cholesky decomposition: A

{tilde over (H)}₀*C⁻¹H₀={tilde over (L)}D{tilde over (L)}* where letting

$A = \begin{bmatrix} e & f^{*} \\ f & h \end{bmatrix}$

we now have

${\overset{\sim}{L} = \begin{bmatrix} 1 & 0 \\ {\overset{\Cup}{l}}_{21} & 1 \end{bmatrix}},{D = \begin{bmatrix} l_{11}^{2} & 0 \\ 0 & l_{22}^{2} \end{bmatrix}},{l_{11}^{2} = e},{{\overset{\Cup}{l}}_{21} = {f/e}},{l_{22}^{2} = {h - {\frac{{f}^{2}}{e}.}}}$

Next, we transform the received signal in (Equation 21) as z₀=D⁻¹{tilde over (L)}⁻¹{tilde over (H)}₀*C⁻¹y₀ and note that z₀ can be expanded as z₀={tilde over (L)}*x₀+{circumflex over (n)}_(o)

$\begin{matrix} {{{E\left\lbrack {{\hat{n}}_{0}{\hat{n}}_{0}^{*}} \right\rbrack} = \begin{bmatrix} {1/l_{11}^{2}} & 0 \\ 0 & {1/l_{22}^{2}} \end{bmatrix}},{{\overset{\sim}{L}}^{*} = \begin{bmatrix} 1 & {\overset{\Cup}{l}}_{21}^{*} \\ 0 & 1 \end{bmatrix}}} & \left( {{Equation}\mspace{14mu} 22} \right) \end{matrix}$

Note that z₀ modeled as in (Equation 22) is precisely the form which allows for square-root-free implementation of two-symbol demodulator derived above (please see Equation 3).

From the above it can be seen that the present invention is advantageous in that it employs a new way of linear transformations and metric computations that are square-root free. The invention employs a look up table based implementation for the metric computations which eliminate redundancy and substantially reduce the number of multiplications. Moreover, inventive method exploits the fact that the un-normalized constellation symbols are complex integers so that the product of a real-number and an un-normalized constellation symbol can be implemented by only additions. The inventive method also enables a greater efficiency for whitening colored noise prior to demodulation, one of which involves no square-root operation. The invention results in less complexity, faster operation, lower power consumption, without sacrificing performance

Having described preferred embodiments of a system and method (which are intended to be illustrative and not limiting), it is noted that modifications and variations can be made by persons skilled in the art in light of the above teachings. It is therefore to be understood that changes may be made in the particular embodiments disclosed which are within the scope of the invention as outlined by the appended claims Having thus described aspects of the invention, with the details and particularity required by the patent laws, what is claimed and desired protected by Letters Patent is set forth in the appended claims. 

1. A method for a square-root free 2 symbol max-log receiver comprising the steps of: i) obtaining linear transformations of a received two stream signal and a channel matrix without implementing square-root operations (100); ii) listing out all possibilities for a first symbol of the said received two stream signal (101); iii) building look-up-tables needed for computing first metrics associated with all possibilities for a first symbol of the said two stream signal (102); iv) determining a second symbol of the two stream signal for each said first symbol listed out (103); v) evaluating said first metrics for each said first symbol and second symbol pair using said look up tables (103); vi) listing out all possibilities for said second symbol 2 (104); vii) building look-up-tables needed for computing second metrics associated with all possibilities for a second symbol of said received two stream signal (105); viii) determining a first symbol for each choice of said second symbol listed out (106); ix) evaluating said second metrics for each said second symbol and first symbol pair using the look up tables (106); vii) determining an exact max-log log likelihood ratio for each coded bit using said second metrics (107); and viii) decoding a least one codeword in the two stream signal using the determined exact max-log log likelihood ratios for all bits (108).
 2. The method of claim 1, wherein said step i) obtaining linear transformations of a received signal comprises obtaining ${\overset{\sim}{Q} = {Q\begin{bmatrix} {1/r_{0,0}} & 0 \\ 0 & {1/r_{1,1}} \end{bmatrix}}},{{\overset{\Cup}{r}}_{0,1} = {r_{0,1}/r_{0,0}}}$ and θ = r_(1, 1)²/r_(0, 0)².
 3. The method of claim 2, wherein said step i) obtaining linear transformations of a received signal comprises multiplying a received vector z with {tilde over (Q)}^(H) giving {tilde over (z)}={tilde over (Q)}^(H)z={tilde over (Q)}^(H)(QRx+n)={tilde over (R)}x+ñ where $\overset{\sim}{z} = {\begin{bmatrix} {\overset{\sim}{z}}^{(0)} \\ {\overset{\sim}{z}}^{(1)} \end{bmatrix} = {{{\begin{bmatrix} 1 & {r_{0,1}/r_{0,0}} \\ 0 & 1 \end{bmatrix}\begin{bmatrix} x^{(0)} \\ x^{(1)} \end{bmatrix}} + \overset{\sim}{n}} = {\begin{bmatrix} {x^{(0)} + {{\overset{\Cup}{r}}_{0,1}x^{(1)}}} \\ x^{(1)} \end{bmatrix} + \begin{bmatrix} {\overset{\sim}{n}}^{(0)} \\ {\overset{\sim}{n}}^{(1)} \end{bmatrix}}}}$ and Q^(H)Q=I with ñ˜N_(c)(0, diag{1r_(0,0) ²,1/r_(1,1) ²}).
 4. The method of claim 1, wherein said step i) obtaining linear transformations of a received signal comprises obtaining {tilde over (Q)},{hacek over (r)}_(0,1), r_(1,1) ² and θ via a modified QR decomposition for any N_(r)≧2 and using a channel matrix {tilde over (H)} and considering {tilde over (H)}=[{tilde over (h)}₀,{tilde over (h)}hd 1] and obtaining α={tilde over (h)}₀ ^(H){tilde over (h)}₁;
 5. The method of claim 1, wherein said step i) obtaining linear transformations of a received signal comprises obtaining {tilde over (Q)},{hacek over (r)}_(0,1), r_(1,1) ² and θ via a modified QR decomposition for any N_(r)≧2 and using a channel matrix {tilde over (H)} and obtaining δ=∥{tilde over (h)}₀∥² and first column of {tilde over (Q)}=[{tilde over (q)}₀,{tilde over (q)}₁] as {tilde over (q)}₀={tilde over (h)}₀/δ.
 6. The method of claim 1, wherein said step i) obtaining linear transformations of a received signal comprises obtaining {tilde over (Q)},{hacek over (r)}_(0,1), r_(1,1) ² and θ via a modified QR decomposition for any N_(r)≧2 and using a channel matrix {tilde over (H)} and obtain the vector {tilde over (v)}={tilde over (h)}₁−α{tilde over (q)}₀; and r_(1,1) ²=∥{tilde over (v)}∥² and the second column of {tilde over (Q)}=[{tilde over (q)}₀,{tilde over (q)}₁] as {tilde over (q)}₁={tilde over (v)}/r_(1,1) ².
 7. The method of claim 1, wherein said step i) obtaining linear transformations of a received signal comprises obtaining {tilde over (Q)},{hacek over (r)}_(0,1), r_(1,1) ² and θ via a modified QR decomposition for any N_(r)≧2 and using a channel matrix {tilde over (H)} and obtaining {hacek over (r)}_(0,1) as {hacek over (r)}_(0,1)=α/δ.
 8. The method of claim 1, wherein said step i) obtaining linear transformations of a received signal comprises obtaining {tilde over (Q)},{hacek over (r)}_(0,1), r_(1,1) ² and θ via a modified QR decomposition for any N_(r)2 and using a channel matrix {tilde over (H)} and obtaining θ as θ=r_(1,1) ²/δ.
 9. The method of claim 1, wherein said steps iv), v) or viii), ix) comprises generating a soft bit k:k=0, 1, . . . B^((v))−1, corresponding to transmitted symbol x^((v)): v=0, 1 using max-log approximation of LLR as follows: ${L\left( b_{k,v} \right)} \approx {{\min\limits_{c_{j} \in S_{k,v,0}}{{\left( {\overset{\sim}{z} - {\overset{\sim}{R}c_{j}}} \right)^{H}\begin{bmatrix} r_{0,0}^{2} & 0 \\ 0 & r_{1,1}^{2} \end{bmatrix}}\left( {\overset{\sim}{z} - {\overset{\sim}{R}c_{j}}} \right)}} - {\min\limits_{c_{j} \in S_{k,v,1}}{{\left( {\overset{\sim}{z} - {\overset{\sim}{R}c_{j}}} \right)^{H}\begin{bmatrix} r_{0,0}^{2} & 0 \\ 0 & r_{1,1}^{2} \end{bmatrix}}\left( {\overset{\sim}{z} - {\overset{\sim}{R}c_{j}}} \right)}}}$
 10. The method of claim 1, wherein said steps iv), v) or viii), ix) comprises generating a soft bit k:k=0, 1, . . . B^((v))−1, corresponding to transmitted symbol x^((v)): v=0, 1 using max-log approximation of LLR wherein a vector c_(j)=[c_(j) ⁽⁰⁾, c_(j) ⁽¹⁾]^(T) with each element c_(j) ^((v)) is a complex symbol in the modulation constellation applied to layer v, and S_(k,v,1) and S_(k,v,0) denote the set containing all possible vector c_(j) such that bit k of element c_(j) ^((v)) is equal to 1 and 0 respectively.
 11. The method of claim 1, wherein said steps iv), v) or viii), ix) comprises considering soft bit calculation for a second layer (i.e. v=1), performing for each c_(j) ⁽¹⁾εS⁽¹⁾:j=0, 1, . . . , K⁽¹⁾−1, computing d_(j) ⁽¹⁾=θ|{tilde over (z)}⁽¹⁾−c_(j) ⁽¹⁾|².
 12. The method of claim 1, wherein said steps iv), v) or viii), ix) comprises considering soft bit calculation for a second layer (i.e. v=1), performing for each c_(j) ⁽¹⁾εS⁽¹⁾: j=0, 1, . . . , K⁽¹⁾−1, computing soft estimate {tilde over (z)}⁽⁰⁾−{hacek over (r)}_(0,1)c_(j) ⁽¹⁾ and then selecting a symbol ĉ_(j) ⁽⁰⁾εS⁽⁰⁾ that is closest to {tilde over (z)}⁽⁰⁾−{hacek over (r)}_(0,1)c_(j) ⁽¹⁾.
 13. The method of claim 1, wherein said steps iv), v) or viii), ix) comprises considering soft bit calculation for a second layer (i.e. v=1), performing for each c_(j) ⁽¹⁾εS⁽¹⁾:j=0, 1, . . . , K⁽¹⁾−1, computing d_(j) ⁽⁰⁾=|{tilde over (z)}⁽⁰⁾−ĉ_(j) ⁽⁰⁾−{hacek over (r)}_(0,1)c_(j) ⁽¹⁾|².
 14. The method of claim 1, wherein said steps iv), v) or viii), ix) comprises considering soft bit calculation for a second layer (i.e. v=1), performing for each c_(j) ⁽¹⁾εS⁽¹⁾:j=0, 1, . . . , K⁽¹⁾−1, computing total distance d_(j) ^(TOTAL)=d_(j) ⁽⁰⁾+d_(j) ⁽¹⁾.
 15. The method of claim 1, wherein said steps iv), v) or viii), ix) comprises: a) Determining ${{\hat{d}}^{(1)} = {\min\limits_{c_{j}^{(1)} \in S^{(1)}}\left\{ d_{j}^{TOTAL} \right\}}},{{\hat{c}}_{j}^{(1)} = {\underset{c_{j}^{(1)} \in S^{(1)}}{\arg \; \min}\left\{ d_{j}^{TOTAL} \right\}}}$ and noting that ĉ_(j) ⁽¹⁾εS⁽¹⁾ is the symbol of layer-1 in the maximum likelihood (ML) decision. Letting {circumflex over (b)}_(k,1):k=0, 1, . . . , B⁽¹⁾−1 be the bit-label of ĉ_(j) ⁽¹⁾εS⁽¹⁾; b) computing each soft bit L(b_(k,1)):k=0, 1, . . . , B⁽¹⁾−1 as ${L\left( b_{k,1} \right)} = {{{r_{0,0}^{2}\left( {{\hat{d}}^{(1)} - {\min\limits_{d_{m}^{TOTAL} \in D_{k}^{1}}d_{m}^{TOTAL}}} \right)}\mspace{14mu} {if}\mspace{14mu} {\hat{b}}_{k,1}} = 0}$ or ${{L\left( b_{k,1} \right)} = {{{r_{0,0}^{2}\left( {{\min\limits_{d_{m}^{TOTAL} \in D_{k}^{0}}d_{m}^{TOTAL}} - {\hat{d}}^{(1)}} \right)}\mspace{14mu} {if}\mspace{14mu} {\hat{b}}_{k,1}} = 1}};$ and c) where D_(k) ¹ and D_(k) ⁰ denote the set containing all distances d_(j) ^(TOTAL) correspond to c_(j) ⁽¹⁾ that has that bit k equal to 1 and 0, respectively.
 16. The method of claim 1, wherein said steps iii) or vii) comprises looking at a max-log LLR expression and after expanding d_(j) ⁽⁰⁾, d_(j) ⁽¹⁾, we drop the terms that do not influence the LLR.
 17. The method of claim 1, wherein said steps iii) or vii) comprises pulling out normalizing factors so that the modulation symbols are from un-normalized QAM constellations which contain integers. Then we can use the fact that the product of an integer and a real number can be obtained using just additions.
 18. The method of claim 1, wherein said steps iii) or vii) comprises forming LUTs containing common terms that mostly depend only on the channel coefficients. The entries of these LUTs are repeatedly accessed while computing all d_(j) ⁽⁰⁾, d_(j) ⁽¹⁾ and multiplications are avoided.
 19. The method of claim 1, wherein said steps iii) or vii) comprises computing d_(j) ⁽¹⁾=θ(c_(j,R) ⁽¹⁾)²+θ(c_(j,I) ⁽¹⁾)²−2{tilde over (z)}_(R) ⁽¹⁾c_(j,R) ⁽¹⁾θ−2{tilde over (z)}_(I) ⁽¹⁾c_(j,I) ⁽¹⁾θ.
 20. The method of claim 1, wherein said steps iii) or vii) comprises computing d_(j) ⁽¹⁾=θ|{tilde over (z)}⁽¹⁾−c_(j) ⁽¹⁾|² in a manner where all {d₀ ⁽¹⁾} can be computed with only 20 real multiplications. 